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Abstract 

We study water flow computation on imprecise terrains. We consider two approaches to 
modeling flow on a terrain: one where water flows across the surface of a polyhedral terrain 

psj in the direction of steepest descent, and one where water only flows along the edges of a 

predefined graph, for example a grid or a triangulation. In both cases each vertex has an 
imprecise elevation, given by an interval of possible values, while its {x, y)-coordinates are 

^ ^ fixed. For the first model, we show that the problem of deciding whether one vertex may be 

contained in the watershed of another is NP-hard. In contrast, for the second model we give 
a simple O(nlogn) time algorithm to compute the minimal and the maximal watershed of 
a vertex, or a set of vertices, where n is the number of edges of the graph. On a grid model, 
we can compute the same in 0{n) time. 

^ Rose knew almost everything that water can do, 

there are an awful lot when you think what. 

Gertrude Stein, The World is Round. 

1 Introduction 

(N 

Simulating the flow of water on a terrain is a problem that has been studied for a long time in 
geographic information science (gis), and has received considerable attention from the compu- 
tational geometry community due to the underlying geometric problems [T| I20| 17]. It can be an 
important tool in analyzing flash floods for risk management [2], for stream flow forecasting [17J, 
^ and in the general study of geomorphological processes [5], and it could contribute to obtaining 

more reliable climate change predictions |26| . 
T-H When modeling the flow of water across a terrain, it is generally assumed that water flows 

downward in the direction of steepest descent. It is common practice to compute drainage 
networks and catchment areas directly from a digital elevation model of the terrain. Most hy- 
drological research in GIS models the terrain surface with a grid in which each cell can drain to 
one or more of its eight neighbors (e.g. [25]). This can also be modeled as a computation on 
a graph, in which each node represents a grid cell and each edge represents the adjacency of 
two neighbors in the grid. Alternatively, one could use an irregular network in which each node 
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Figure 1: Left: An imprecise terrain. Each vertex of the triangulation has a elevation interval 
(gray). Center: a realization of the imprecise terrain. Right: the same realization together with 
the highest and lowest possible realizations of the imprecise terrain. 

drains to one or more of its neighbors, which may reduce the required storage space by allowing 
less interesting parts of the terrain to have a lower sampling density. We will refer to this as the 
network model, and we assume that, from every node, water flows down along the steepest inci- 
dent edge. Assuming the elevation data is exact, drainage networks can be computed efficiently 
in this model (e.g. |6j]). In computational geometry and topology, researchers have studied flow 
path and drainage network computations on triangulated polyhedral surfaces (e.g. [8l[9l[19]). In 
this model, which we call the surface model, the flow of water can be traced across the surface 
of a triangle. This avoids creating certain artifacts that arise when working with grid models. 
However, the computations on polyhedral surfaces may be significantly more difficult than on 
network models [S]. 

Naturally, all computations based on terrain data are subject to various sources of uncer- 
tainty, like measurement, interpolation, and numerical errors. The GIS community has long 
recognized the importance of dealing with uncertainty explicitly, in particular for hydrologi- 
cal modeling. A common approach is to model the elevation at a point of the terrain using 
stochastic methods [28]. However, the models available in the hydrology literature are unsatis- 
factory [3l [231 [21] and computationally expensive [27]. A particular challenge is posed by the 
fact that hydrological computations can be extremely sensitive to small elevation errors [14^ [T8] . 
While most of these studies have been done in the network model, we note that there also exists 
work on the behaviour of watersheds on noisy terrains in the surface model by Haverkort and 
Tsirogiannis |13j . 

A non-probabilistic model of imprecision that is often used in computational geometry con- 
sists in representing an imprecise attribute (such as location) by a region that is guaranteed to 
contain the true value. This approach has also been applied to polyhedral terrains (e.g. [12^116]). 
replacing the exact elevation of each surface point by an imprecision interval (see Figure [l]). In 
this way, each terrain vertex does not have one fixed elevation, but a whole range of possible el- 
evations which includes the true elevation. Choosing a concrete elevation for each vertex results 
in a realization of the imprecise terrain. The realization is a (precise) polyhedral terrain. Since 
the set of all possible realizations is guaranteed to include the true (unknown) terrain, one can 
now obtain bounds on parameters of the true terrain by computing the best- and worst-case 
values of these parameters over the set of all possible realizations. Note that we assume the 
error only in the z-coordinate (and not in the x, y-coordinates). This is partially motivated by 
the fact that commercial terrain data suppliers often only report elevation error |1U| . However, 
it is also a natural simplification of the model, since the true terrain needs to have an elevation 
at any exact position in the plane. 

In this paper we apply this model of imprecise terrains to problems related to the simulation 
of water flow, both on terrains represented by surface models and on terrains represented by 
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network models. One of the most fundamental questions one may ask about water flow on 
terrains is whether water flows from a point p to another point q. In the context of imprecise 
terrains, reasonable answers may be "definitely not", "possibly", and "definitely". The water- 
shed of a point in a terrain is the part of the terrain that drains to this point. Phrasing the 
same question in terms of watersheds leads us to introduce the concepts of potential (maximal) 
and persistent (minimal) watersheds. 



Results In Section [3] we show that the problem of deciding whether water can flow between 
two given points in the surface model is NP-hard. Fortunately, the situation is much better for 
the network model, and therefore as a special case also for the D-8 grid model which is widely 
adopted in GIS applications. In Section [4] and Section [5] we present various results using this 
model. In Section [4. 1| we present an algorithm to compute the potential watershed of a point. 
On a terrain with n edges, our algorithm runs in O(nlogn) time; for grid models the running 
time can even be improved to 0{n). We extend these techniques and achieve the same running 
times for computing the potential downstream area of a point in Section 4.2 and its persistent 
watershed in Section 14.31 In order to be able to extend these results in the network model, 



we define a certain class of imprecise terrains which we call regular in Section [5j We prove 



that persistent watersheds satisfy certain nesting conditions on regular terrains in Section 5.3 



This leads to efficient computations of fuzzy watershed boundaries in Section 5.4 and to the 
definition of the fuzzy ridge in Section 5.5, which delineates the persistent watersheds of the 
"main" minima of a regular terrain and which is equal to the union of the areas where the 
potential watersheds of these minima overlap. We can compute this structure in 0(n log n) 
time (see Theorem [g]) and we discuss an algorithm that turns a non-regular terrain into a 
regular one (see Section 5.2). We conclude the paper with the discussion of open problems in 
Section [H 



2 Preliminaries 

We give the definition of imprecise terrains and realizations and discuss the two fiow models 
used in this paper. 

2.1 Basic definitions and notation 

We define an imprecise terrain T as a possibly non-planar geometric graph G with nodes 
V C and edges C y x where each node v ^ V has an imprecise third coordinate, 
which represents its elevation. We denote the bounds of the elevation of v with low{v) and 
high{v). A realization R of an imprecise terrain T consists of the given graph together with 
an assignment of elevations to nodes, such that for each node v its elevation elevji{v) is at least 
low{v) and at most high{v). As such, it is a fully embedded graph R = {Vr^Eji) in IR'^, where 
Vr is defined as the set {(x,y,z) | 3 u G 1/ : v = (x, y), z = elevii{v)}. Note that this defines a 
one-to-one correspondence between nodes of V and Vr. The edge set Er is induced by E under 
this correspondence. With slight abuse of notation we will sometimes refer to nodes of V by their 
corresponding nodes in Vr. We denote with R~ the realization, such that elevR-{v) = low{v) 
for every vertex v and similarly the realization it!"'", such that elevR+{v) = high{v). The set of 
all realizations of an imprecise terrain T is denoted with 'Rt- 

Now, consider a realization R of an imprecise terrain as defined above. For any set of nodes 
P Q Vr, we define the neighborhood of P as the set N{P) = {s : s ^ P A3 t £ P : {s,t) £ Er}. 
If P is a connected set, all nodes of P have the same elevation and this elevation is strictly lower 
than the elevation of any node in N{P), then P constitutes a local minimum. Likewise, a 



3 



local maximum is a set of nodes at the same elevation of which the neighborhood is strictly 
lower. 

2.2 A model of discrete water flow 

Consider a realization R of an imprecise terrain as defined above. If water is only allowed to 
flow along the edges of the realization, then the realization represents a network. Therefore 
we refer to this model of water flow as the network model. Below, we state more precisely 
how water flows in this model and give a proper definition of the watershed. This model or 
variations of it have been used before, for example in [6l [22l I25| . 

The steepness of descent (slope) of an edge {p,q) £ Er is defined as crji{p,q) = {elevB,{p) — 
devji(q))/\pq\, where \pq\ is the Euclidean distance between the corresponding nodes in V. The 
node g is a steepest descent neighbor of p, if and only if (tr{p, q) is non-negative and maximal 
over all neighbours q of p. Water that arrives in p will continue to flow to each of its steepest 
descent neighbors, unless p constitutes a local minimum. If there exists a local minimum P 3 p, 
then the water that arrives in p will flow to the neighbors of p in P and eventually reach all the 
nodes of P, but it will not flow further to any node outside the set P. If water from p reaches 
a node g G Vr then we write Pj^q ( "p flows to q in R" ) , and for technical reasons we define 
Pj^p for all p and R. 

The discrete watershed of a node g in a realization R is defined as the union of nodes 
that flow to q in R, that is Wniq) ■= {p ■ Pjf^q}- Similarly, we define the discrete watershed 
of a set of nodes Q in this realization as W/{(Q) := Ugeg'^nC'?)- 

Consider the graph G of the imprecise terrain. Let vr be a path in G, with no repeated 
vertices. We say vr is a flow path in a realization R if it carries water to a local minimum 
and visits all nodes of this local minimum in R. For two consecutive vertices p, g in vr, g is a 
steepest descent neighbour of p. For any pair of nodes p, q in vr, we write p j^q if pj^ q, that is, 
vr contains p and q in this order. We denote with vr[p, q] the subpath of vr from p to q, including 
these two nodes. For any given set of realizations S C 3^7^, we denote with n(S') the set of all 
flow paths in any realization in S. 

2.2.1 Flow paths are stable 

This subsection is a note on flow paths, which we defined for the network model above. We 
define when a flow path is stable and argue that any flow path induced by a realization in JIt 
is stable with respect to some e— neighborhood of JIt- Intuitively, the analysis in this section 
shows that the flow paths considered in our model are never the result of an isolated degenerate 
situation, but could also exist if the estimated elevation intervals of the vertices would be slightly 
different. This may serve as a justification or proof of soundness of the network model. 

If for two realizations R, R' G and any node v £ Vr and its corresponding node v' G Vri 
we have \vv'\ < e, then we call R! an e— perturbation of R. For a set of realizations 5", let 

denote the union of S with the e— perturbations of elements of S. We say that a flow 
path vr G n(S') is stable with respect to S if for some e > the flow path exists in any 
e— perturbation of some R £ S (we call R the perturbation center). Let 11(5) C 11(5') denote 
the subset of flow paths that is stable with respect to S. We call a realization which does 
not contain horizontal edges and in which any node has at most one steepest descent neighbor 
non- ambiguous , similarly, a realization for which any of these properties does not hold is 
called ambiguous. We have the following lemma, which implies that any flow path induced by 
an element of JIt is stable with respect to Jlj^, for any e > 0. 

Lemma 1 For any set of realizations S, we have that 11(5) C lime_j.o 11(5^). 
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Proof: Given any value e > 0, we argue that the set 11(5) is contained in the set 11(5'^). 
Clearly, any flow path induced by a non-ambiguous realization R is stable with perturbation 
center R. Now, let vr = pi,p2, ■ ■ ■ ,Pk be a flow path from pi to pk which is induced by an 
ambiguous realization R £ S. We lower each node pi by e/2 + {ie)/{Ak) and perturb the 
remaining vertices by some value smaller than e/4 to create a non-ambiguous realization which 
also induces vr. This proves the claim. □ 



2.3 A model of continuous water flow 

Consider an imprecise terrain, where the graph that represents the terrain forms a planar 
triangulation in the {x, y)-domain. Any realization of this terrain is a polyhedral terrain with a 
triangulated surface. If we assume that the water which arrives at a particular point p on this 
surface will always flow in the true direction of steepest descent at p across the surface, possibly 
across the interior of a triangle, then we obtain a continuous model of water flow. Since the 
steepest descent paths do not necessarily follow along the edges of the graph, but instead lead 
across the surface formed by the graph, we call this model the surface model. This model has 
been used before, for example in [HI [9l [T9]. 

Since, as we will show in the next section, it is already NP-hard to decide whether water 
from a point p can potentially flow to another point q, we will focus on the network model in 
the rest of the paper, and we do not formally define watersheds in the surface model. Therefore, 
we will simply use the term "watershed" to refer to discrete watersheds in this paper. 



3 NP-hardness in the surface model 



In the surface model water flows across the surface of a polyhedral terrain; refer to Section 2.3 



for the details of the model. In this section we prove that it is NP-hard to decide whether water 
potentially flows from a point s to another point t in this model. The reduction is from 3-SAT; 
the input is a 3-CNF formula with n variables and m clauses. We first present the general idea 
of the proof, then we proceed with a detailled description of the construction, and finally we 
prove the correctness. 



3.1 Overview of the construction 

The main idea of the NP-hardness construction is to encode the variables and clauses of the 
3-SAT instance in an imprecise terrain, such that a truth assignment to the variables corresponds 
to a realization — i.e., an assignment of elevations — of this terrain. If and only if all clauses are 
satisfied, water will fiow from a certain starting vertex s to a certain target vertex t. We first 
introduce the basic elements of the construction: channels and switch gadgets. 



Channels We can mold channels in the fixed part of the terrain 
to route water along any path, as long as the path is monotone 
in the direction of steepest descent on the terrain. We do this 
by increasing the elevations of vertices next to the path, thus 
building walls that force the water to stay in the channel. We 
can end a channel in a local minimum anywhere on the fixed part 
of the terrain, if needed. 




Switch gadgets The general idea of a switch gadget is that it provides a way for water to 
switch between channels. A simple switch gadget has one incoming channel, three outgoing 
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channels, and two control vertices a and b, placed on the boundary of the switch. The water 
from the incoming channel has to flow across a central triangle, which is connected to a and b. 
Depending on their elevations, the two vertices a and b divert the water from the incoming 
channel to a particular outgoing channel and thereby "control" the behaviour of the switch 
gadget. This is possible, since the slope of the central triangle, which the water needs to 
pass, depends on the elevations of a and b and those two are the only vertices with imprecise 
elevations. The elevations of the vertices which are part of the channels are fixed. Refer to 
Figure [2] for an illustration. 




Figure 2: Three different states of a simplistic switch gadget. 



We can also build switches for multiple incoming channels. In this case, every incoming 
channel has its own dedicated set of outgoing channels, and it is also controlled by only two 
vertices, see Figure [3} Note that we can lead the middle outgoing channel to a local minimum 
as shown in the examples and, in this way, ensure that, if any water can pass the switch, the 
elevations of its control vertices are at unambigous extremal elevations. Depending on the 
exact construction of the switch, we may want them to be at opposite extremal elevations or at 
corresponding extremal elevations. 




Figure 3: Sketch of a switch with multiple incoming channels. 



Global layout The global layout of the construction is depicted in Figure |4j The construction 
contains a grid of mxn cells, in which each clause corresponds to a column and each variable to 
a row of the grid. The grid is placed on the western slope of a "mountain" ; columns are oriented 
north-south and rows are oriented east- west. We create a system of channels that spirals around 
this mountain, starting from s at the top and ending in t at the bottom of the mountain. We 
ensure that in no realization, water from s can escape this channel system and, if it reaches t, 
we know that it followed a strict course that passes through every cell of the grid exactly once, 
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Figure 4: Left: Global view of the NP-hardness construction, showing the grid on the mountain 
slope. The fixed parts are shown in gray, the variable parts are shown light yellow and the 
imprecise vertices are filled light green; Right: Detail of a clause, which forms one of the 
columns of the grid. 

column by column from east to west, and within each column, from north to south. Embedded 
in this channel system, we place a switch gadget in every cell of the grid, which allows the 
water from s to "switch" from one channel to another within the current column depending on 
the elevations of the vertices that control the gadget. In this way, the switch gadgets of a row 
encode the state of a variable. To ensure that the state of a variable is encoded consistently 
across a row of the grid, the switch gadgets in a row are linked by their control vertices. Every 
column has a dedicated entry point at its north end, and a dedicated exit point at its south 
end. If and only if water flows between these two points, the clause that is encoded in this 
column is satisfied by the corresponding truth assignment to the variables. The slope of the 
mountain is such that columns descend towards the south, and the exit point of each column 
(except the westernmost one) is lower than the entry point of the adjacent column to the west; 
water can flow between these points through a channel around the back of the mountain. The 
easternmost column's entry point is the starting vertex s, and the westernmost column's exit 
point is the target vertex t. 

Clause columns To encode each clause, we connect the switch gadgets in a column of the 
global grid by channels in a tree-like manner. By construction, water will arrive in a different 
channel at the bottom of the column for each of the eight possible combinations of truth values 
for the variables in the clause. This is possible because a switch gadget can switch multiple 
channels simultaneously. We let the channel in which water would end up if the clause is not 
satisfied lead to a local minimum; the other seven channels merge into one channel that leads 
to the exit point of the clause. The possible courses that water can take will also cross switch 
gadgets of variables that are not part of the clause: in that case, each course splits into two 
courses, which are merged again immediately after emerging from the switch gadget. Figure |4] 
(right) shows an example. 

Sloped switch gadgets Since the grid is placed on the western slope of a mountain, water 
on the central triangle of a switch will veer off towards the west, regardless of the elevations of 
its control vertices. However, as we will see, we can still design a working switch gadget in this 
case. Recall that we link the switch gadgets of a variable row by their control vertices, such 
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that each switch gadget shares one control vertex with its neighboring cell to the west and one 
with its neighboring cell to the east. As mentioned before, such a row encodes the state of a 
particular variable. We say that it is in a consistent state if either all control vertices of the 
switches are high or all control vertices are low. Thus, we will use the following assignment of 
truth values to the elevations of the control vertices of our switch gadgets: both vertices set to 
their highest elevation encodes true; both vertices set to their lowest elevation encodes false; 
other combinations encode confused. Depending on the truth value encoded by the elevations of 
the imprecise vertices, water that enters the gadget will flow to different channels. The channels 
in which the water ends up when the gadget reads confused always lead to a local minimum. 
For the other channels, their destination depends on the clause. In Figure[5]you can see a sketch 
of a sloped switch gadget which works the way described above. 




Figure 5: Illustration of a sloped switch gadget similar to the one used in the final construction. 
The final gadget has multiple incoming channels, which is not shown in this figure. 

3.2 Details of the construction 

Recall that we are given a 3-SAT instance with n variables and m clauses. The central part 
of the construction, which will contain the gadgets, consists of a grid of n rows — one for each 
variable — and m columns — one for each clause. We denote the width of each row, measured 
from north to south, hj B = 400, and the width of a column, measured from west to east, by 
A = max((n + 1) B, 4000). Ignoring local details, on any line from north to south in this part of 
the construction, the terrain descends at a rate of dz/dy = 1, and on any line from east to west, 
it descends at a rate of dz/dx = 1; thus we have z = x + y. Observe that each column measures 
nB < A from north to south; thus the southern edge of each column is at a higher elevation 
than the northern edge of the next column to the west. The dedicated entrance and exit points 
of column 1 < j < n are placed at {jA — ^A,nB,jA — ^A + nB) and {jA — ^A,0,jA — ^A), 
thus allowing the construction of a descending channel from each column's exit point to the 
entry point of the column to the west. 

For every variable Vi, 1 < i < n, we place m + 1 imprecise vertices Vij, for < j < n, 
in row i, on the boundaries of the columns corresponding to the m clauses. Vertex Vij has x- 
coordinate jA, y-coordinate iB — ^B, and an imprecise z-coordinate [jA + iB — ^B, 



8 




Figure 6: Distances and gradients on a connector gadget. All coordinates are relative to the 
lowermost position of the control vertex d = f The other control vertex is e = Vij. Thus, 

e and d are the only imprecise vertices. The x- and y-coordinates of the vertices are indicated 
on the axes. The elevations of the key vertices are written next to the vertices. The elevations 
of the control vertices are expressed as a function of G [0, 1]. The directions of steepest 
descent on the different faces of the gadget are expressed in the form dy/dx, as a function of a 
and /3. 

jA + iB — \B + 20]. On every pair of imprecise vertices Vi(^j_i^,Vij we build a switch gad- 
get Gij] thus there is a switch gadget for each variable/clause pair. The coordinates of the 
vertices in each gadget, relative to the coordinates of can be found in Figure |6j 

Switch gadget construction We use the sloped switch gadget described above and illus- 
trated in Figure [5j Our switch gadget occupies a rectangular area that is A wide from west to 
east, and 41 wide from north to south. Its key vertices and their coordinates, relative to each 
other, can be found in Figure [6| There are two imprecise vertices, d and e, with elevation range 
[0,20] and [A, A + 20], respectively — so in any realization, their elevations have the form 20a 
and A + 20/3, respectively, where a, (3 £ [0, 1]. 

On the north edge of the gadget, there may be many more vertices, all collinear with a, b 
and c. The vertices on the western half of the north edge are connected to the western control 
vertex, and the vertices on the eastern half of the north edge are connected to the eastern 
control vertex. In particular, each gadget Gij is designed to receive water from four channels 
that arrive at four points Siji, Sij2, Sij3, Sij^ on the north edge, close to b; the coordinates of 
these points are stjk = {^A - 150 + GOk, 40, ^A - 110 + 60k). 

On the south edge of the gadget, there is a similar row of vertices, all collinear with /, 
g and h, that are connected to the control vertices. To the south, the gadget is connected 
to twelve channels that catch all water that arrives at certain intervals on the south edge: 
for each k G {1,2,3,4}, there is a western channel tijk catching all water arriving between 
Sijk — (82,41,123) and sijk — (77,41,118), a middle channel Cijk catching all water arriving 
between Sijk — (77, 41, 118) and s-ijk — (44, 41, 85), and an eastern channel fijk catching all water 
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arriving between Sjj/j — (44,41,85) and Sijk - (39,41,80). 

In a particular realization R, we define the switch gadget to be in a false state if q = /? = 0, 
in a true state if a = /3 = 1, and in a confused state if a < ^ while (3 > ^, or if a > ^ while 
/3 < |. As we will show below, in the true, false, and confused states the gadget leads any water 
that comes in at any point Sijk into tijk, fijk, and Cijk, respectively. 

We model the fixed part of the terrain such that the middle channels all lead to local minima. 
The western and eastern channels correspond to a (partial) truth assignment of the variables of 
the clause that is represented by the column that contains the gadget; these channels lead to a 
local minimum or to the next row, as described below. 

Constructing the clause columns Each clause is modeled in a column j by making certain 
connections between the outgoing channels of each gadget to the dedicated entrance points of 
the gadget in the next row. Observe that by our choice of B, the entrance point of column j 
lies above all entrance points of Gnj, all outgoing channels of any gadget G(j_|_i)j start at higher 
elevations than all entrance points of Gij, and all outgoing channels of Gij start at an elevation 
higher than the exit point of the column. This ensures that all channels described below can 
indeed be built as monotonously descending channels, so that water can flow through it. We 
will now explain the connections which we use to build a clause. 

Let p > q > r he the indices of the variables that appear in the clause. The water courses 
modelling the clause start at the entry point of the column, from which any water is led through 
a channel to entry point Snji of gadget Gnj- 

For i / {p, q, r}, k £ {1, 2, 3, 4}, we connect both tijk and fijk to S(j_i)jfc (if i > 1) or to the 
exit point of the column (if i = 1). 

We connect tpji and fpji to S(p_i)ji and S{p_i)j2) respectively. Thus, for i S {q, ...,p — 1}, 
water that enters Gij at Siji and Sij2 represents the cases that p is true and p is false, respectively. 

We connect tgji, fqji,tgj2 and fgj2 to S(g_i)ji, S(g_i)j2, S(g-i)j3 and S(g_,i)j4, respectively. 
Thus, for i G {r,...,q — 1}, water that enters Gij at Siji, Sij2, Sijs and Sij^ represents the four 
different possible combinations of truth assignments to p and q, respectively. 

The eight channels tj-ji, frji,trj2i frj2,trj3, frj3,trji, frj4: now represent the eight different 
possible combinations of truth assignments to the variables of the clause. The channel that 
corresponds to the truth assignment that renders the clause false, is constructed such that it 
ends in a local minimum. The other seven channels all lead to (if r > 1) or to the exit 

point of the clause column (if r = 1). 

Analysis of flow through a gadget Below we will analyse where water may leave a gadget 
Gij after entering the gadget at point Sijk, with x-coordinate x^. In the discussion below, all 
coordinates are relative to the lowermost position of the western control vertex of the gadget — 
refer to Figure [6| which also shows the directions of steepest descent (expressed as dx/dy) on 
each face of the gadget. 

First observe that in any case, the directions of steepest descent on Aabd, Abce and Abed 
are at least 1 - 20/A > 199/200 = 0.995 and at most (1 + 20/A)/(l/2) < 201/100 = 2.01. 
Thus, when the water reaches y-coordinate 38, it will be at x-coordinate at least x^ — 4.02 and 
at most Xfc — 1.99. 

Note that the line bd intersects the plane y = 38 at x = — < ^A — 100, and the line 
be intersects the plane ?/ = 38atx = ^A + -^A > ^^ + 100. By our choice of coordinates for the 
entrance points Sijk, we have |xfc — < 90; therefore the water will be on Abed when it reaches 
y = 38. Let ^max and gmm be the maximum and minimum possible gradients dx/dy on Abed, 
respectively. Thus, the water will reach the line de at x-coordinate at least x^ — 4.02 — 385(max 
and at most x^ — 1.99 — 38ymin- 
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Finally, the directions of steepest deseent on Adgf, Aegh and Adeg are more than and 
less than 1 + 20/^1 < 201/200 < 1.01. Thus, the water will reach the line fh at x-coordinate 
more than Xk — 5.03 — SS^max and less than Xk — 1.99 — 38gmin- 

We will now consider five classes of configurations of the control vertices in the gadget, and 
compute the interval of x-coordinates where water may reach the line fh in each case. 

• a = (3 = 1 (true state) In this case we have ^max = ffmin = 2, so water will reach the 
line fh within the x-coordinate interval (xfe — 81.03, x^ — 77.99), and thus it will flow into 
channel tij^- 

• a + P > I (true-ish state) In this case we have gimax < (1 + 20/A)/(l/2) < 2.01 and 
5min > (1 - 20/^)/(l - 3/8) > 199/125 > 1.59. Thus water will reach the line fh within 
the x-coordinate interval (xfe — 81.41, Xfe — 62.41), and thus it will flow into channel tijk 
or Cijk- 

• |<a + /?<| (this includes all confused states) In this case we have g'max > (1 + 
20/A)/(l - 3/8) < 201/125 < 1.61 and gmin > (1 - 20/A)/(l - 1/8) > 199/175 > 1.13. 
Thus water will reach the line fh within the x-coordinate interval (x^ — 66.21, Xfe — 44.93), 
and thus it will flow into channel Cijk- 

• a + P < \ (false-ish state) In this case we have ^max < (1 + 20/A)/(l - 1/8) < 201/175 < 
1.15 and ^min > (1 - 20/A) > 199/200 > 0.99. Thus water will reach the hne fh within 
the x-coordinate interval (x^ — 48.73, x^ — 39.61), and thus it will flow into channel Cijk 
or fijk- 

• a = /3 = (false state) In this case we have gmax = ffmin = 1, so water will reach the 
line fh within the x-coordinate interval (xfc — 43.03, Xfe — 39.99), and thus it will flow into 
channel fijk- 

Correctness of the NP-hardness reduction 

Lemma 2 // water flows from s to t in some realization, then there is a truth assignment of 
the variables of the 3- CNF formula that satisfies the formula. 

Proof: Water that starts flowing from s, which is the entrance point of the clause column m, 
is immediately forced into a channel to entrance point Snmi of gadget Gnm- As calculated above, 
any water that enters a gadget at one of its designated entrance points will leave the gadget 
in one of its designated channels, which leads either to a local minimum, or to a designated 
entrance point of the next gadget. Therefore, water from s can only reach t after flowing through 
all switch gadgets. 

Since all middle outgoing channels Cijk lead to local minima, we know that if there is a flow 
path from s to t, then the water from s is nowhere forced into a middle outgoing channel. It 
follows that no gadget is in a confused state. As a consequence, in any row, either all gadgets 
have their control vertices in the lower relatively open half of their elevation range, or all gadgets 
have their control vertices in the upper relatively open half of their elevation range. In the first 
case, all gadgets in the row are in a false-ish state, and any incoming water from s leaves those 
gadgets in the same channels as if the gadgets were in a proper false state. In the second case, 
all gadgets in the row arc in a true-ish state, and any incoming water from s leaves those gadgets 
in the same channels as if the gadgets were in a proper true state. 

We can now construct a truth assignment A to the variables, in which each variable is true 
if the control vertices in the corresponding row are in the upper halves of their elevation ranges. 
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and false otherwise. It follows from the way in which channel networks in clause columns are 
constructed, that in each clause column, water will flow into one of the seven channels that 
corresponds to a truth assignment that satisfies the corresponding clause — otherwise the water 
would not reach t. Therefore, A satisfies each clause, and thus, the complete 3-CNF formula. 

□ 

Lemma 3 If there is a truth assignment to the variables that satisfies the given 3-CNF formula, 
then there is a realization of the imprecise terrain in which water flows from s to t. 

Proof: We set all control vertices in rows corresponding to true variables to their highest 
positions and all control vertices in rows corresponding to false variables to their lowest positions. 
One may now verify that, by construction, in each clause column water from the column's entry 
point will flow into one of the seven channels that lead to the column's exit point, and thus, 
water from s reaches t. □ 



Thus, 3-SAT can be reduced, in polynomial time, to deciding whether there is a realization 
of T such that water can flow from s to t. We conclude that deciding whether there exists a 
realization of T such that water can flow from s to t is NP-hard. 

Theorem 1 Let T be an imprecise triangulated terrain, and let s and t be two points on the 
terrain. Deciding whether there exists a realization R E "Rt such that p^^qis NP-hard. 



4 Watersheds in the network model 

In the network model we assume that water flows only along the edges of a realization. More 
specifically, water that arrives in a node p continues to flow along the steepest descent edges 
incident on p, unless p is a local minimum. For a formal definition of the watershed and flow 



paths please refer to Section 2.2 



4.1 Potential watersheds 

The potential watershed of a set of nodes Q in a terrain T is defined as 

Wu(Q):= U U^«(^)' 

which is the union of the watersheds of Q over all realizations of T. This is the set of nodes for 
which there exists a flow path to a node of Q. With slight abuse of notation, we may also write 
Wu {q) to denote the potential watershed of a single node q. 



4.1.1 Canonical realizations 

We prove that for any given set of nodes Q in an imprecise terrain, there exists a realization 
R such that W/j((5) = Wy (Q). For this we introduce the notion of the overlay of a set of 
watersheds in different realizations of the terrain. Informally, the overlay is a realization that 
sets every node that is contained in one of these watersheds to the lowest elevation it has in any 
of these watersheds. 

Definition 1 Given a sequence of realizations Ri,...,Rk and a sequence of nodes qi,...,qk, the 
watershed-overlay o/W/j^(gi), W/j^. (g^) is the realization R such that for every node v, we 
have that elev-j^{v) = high{v) if v ^ U"^Ri('?«) '^^^ otherwise 

elev^{v) = min elevji^{v). 
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Lemma 4 Let R be the watershed-overlay o/ W/j^ (gi ),..., W/j^ (g^), and let Q = Ui<i<fcftj 
thenW^{Q) contains Wfi-{qi). 

Proof: Let n be a node of the terrain. We prove the lemma by induction on increasing symbohc 
elevation to show that if u is contained in one of the given watersheds, then it is also contained 
in Wj^{Q). We define level{Ri,u) as the smallest number of edges on any path along which 
water flows from u to qi in Ri\ if there is no such path, then level{Ri,u) = oo. Now we define 
the symbolic elevation of u, denoted elev*{u), as follows: if u is contained in any watershed 
Wr. (gj), then elev*{u) is the lexicographically smallest tuple {elevR-{u), level{Ri,u)) over all i 
such that u G Wij.(gi); otherwise elev*{u) = {high{u) , oo) . 

Now consider a node u that is contained in one of the given watersheds. The base case is that 
u is contained in Q, and in this case the claim holds trivially. Otherwise, let Ri be a realization 
such that u G and such that {elevfi.(u), level{Ri,u)) is lexicographically smallest over 

all 1 < i < /c. By construction, we have that elev^-[u) = elev^{u). Consider a neighbour v 
of u such that (u, v) is a steepest-descent edge incident on u in and level{Ri,v) is minimal 
among all such neighbours v of u. Since elevj^{v) < elevji.{v) < elevii^{u) = elev^{u) and 
level{Ri,v) = level{Ri,u) — 1, it holds that v has smaller symbolic elevation than u. Therefore, 
by induction, v G 'W-^(Q). If v is still a steepest descent neighbor of u in R, then this implies 
u £ Wj^{Q). Otherwise, there is a node v such that aj^{u,v) > aj^{u,v) > 0. There must be 
an Rj such that v £ Wptjigj), since otherwise, by construction of the watershed-overlay, we 
have elev^{v) = high(v) > elevR-{v) and thus, aR^{u,v) > a^{u,v) > a^{u,v) > aR^{u,v) 
and V would not be a steepest descent neighbor of u in Ri. Moreover, we have a^{u,v) > 
and, therefore, elev^{v) < elevj^{u), so v has smaller symbolic elevation than u. Therefore, by 
induction, also v G ^]^{Q) and thus, u G Wj^(Q). □ 

The above lemma implies that for any set of nodes Q, the watershed-overlay R of the 
watersheds of the elements of Q in all possible realizations JIt, would realize the potential 
watershed of Q. That is, we have that Wy (Q) C W-^((5) and since Wy (Q) is the union of 
all watersheds of Q in all realizations, we also have that 'W^iQ) ^ (Q), which implies the 
equality of the two sets. Therefore, we call R the canonical realization of the potential 
watershed Wy (Q) and we denote it with Ru{Q). 

Note, however, that it is not immediately clear that the canonical realization always exists: 
the set of possible realizations is a non-discrete set, and thus the elevations in the canonical 
realization are defined as minima over a non-discrete set. Therefore, one may wonder if these 
minima always exist. Below, we will describe an algorithm that can actually compute the 
canonical realization of any set of nodes Q; from this we may conclude that it always exists. 

4.1.2 Outline of the potential watershed algorithm 

Next, we describe how to compute Wy (Q) and its canonical realization RuiQ) for a given set of 
nodes Q. Note that for all nodes p ^ Wy (Q), we have, by definition of the canonical realization, 
dev ji^(^Q^{p) = high(p). The challenge is therefore to compute Wy (Q) and the elevations of the 
nodes of Wy (Q). Below we describe an algorithm that does this. 

The idea of the algorithm is to compute the nodes of Wy (Q) and their elevations in the 
canonical realization in increasing order of elevation, similar to the way in which Dijkstra's 
shortest path algorithm computes distances from the source. The complete algorithm is laid 
out in Algorithm [l} The correctness and running time of the algorithm are proved in Theorem [2] 
A key ingredient of the algorithm is a subroutine, ExPAND(q', z'), which is defined as follows. 

Definition 2 Let Expand (g', z') denote a function that returns for a node q' and an elevation 
z' G [low{q'), high{q')] a set of pairs of nodes and elevations, which includes the pair {p,z) if 
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and only if p €z N{q'), there is a realization R with elevR{q') G [z' , high{q')] such thatp^q', 
and z is the minimum elevation of p over all such realizations R. 



Algorithm 1 PotentialWS((5) 

1: For all q G Q: Enqueue {q, z) with key z = low{q) 

2: while the Queue is not empty do 

3: {q',z') = DequeueMin() 

4: if q' is not already in the output set then 

5: Output q' with elevation z' 

6: Enqueue each (p, z) G ExPAND(g', z') 

7: end if 

8: end while 



4.1.3 Expansion of a node using the slope diagram 

Before presenting the algorithm for the expansion of a node, we discuss a data structure that 
allows us to do this efficiently. 

We define the slope diagram of a node p as the set of points % = {5i, high{qi)), such that 
qi is a neighbor of p and 6i is its distance to p in the (x, y)-projection. Let qi,q2,---, be a 
subset of the neighbors of p indexed such that qi,q2, ... appear in counter-clockwise order along 
the boundary of the convex hull of the slope diagram, starting from the leftmost point and 
continuing to the lowest point. We ignore neighbors that do not lie on this lower left chain. 

Let Hi be the halfplane in the slope diagram that lies above the line through q^ and gi+i. 
Let U{p) be the intersection of these halfplanes Hi,H2, the halfplane right of the vertical line 
through the leftmost point, and the halfplane above the horizontal line through the bottommost 
point of the convex chain, see the shaded area in Figure [7j Let Zi be the z-value of the point 
where the line through % and %+! intersects the vertical axis of the slope diagram. 

Note that, out of the neighbors of p, each set to their highest position, qi is a steepest descent 
neighbor if and only if the elevation of p lies in the interval {zi, Zi-i). This observation will help 
us to compute the minimum elevation of p such that water flows to any particular neighbor in 
0{\ogd) time, given that we have the slope diagram of p at hand. 

For a neighbor p of q' , we can now compute the elevation of p as it should be returned 
by ExPAND(g', z') by computing the lower tangent to U{p) which passes through the point 
q' = {6',z'), where 6' is the distance from q' to p in the (x, ?/)-projection. This can be done 



high{qi) 
Zi-l Y> 



high{qi) 
k 






Figure 7: Left: slope diagram. Right: querying the slope diagram. 
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via a binary search on the boundary of U{p). Intuitively, this tangent intersects the corner of 
U{p) which corresponds to the neighbor of p that the node q' has to compete with for being the 
steepest-descent neighbor of p. The elevation z at which the tangent intersects the vertical axis, 
is the lowest elevation of p such that q' does not lose, see the figure. In proving the following 
lemma we describe the details of this procedure more specifically. 

Lemma 5 Given the slope diagrams of the neighbours of q' , we can compute the function 
Expand (g', z') in time 0{dlogd') , where d is the node degree of q' , d' is the maximum node 
degree of a neighbor of q' , and n is the number of edges of the terrain. 

Proof: Let p he a neighbor of q' and let Zmin be max{low{p), z'). Obviously, Zmin is a lower 
bound on the elevation that p could have while still allowing flow to q' . There are three cases 
for the outcome of a query with q' in the slope diagram of p. 

(i) If q' lies in the interior of U{p), then q' can never be a steepest descent neighbor of p in a 
non-ambiguous realization. As such, p is not included in the result of Expand (g', z'). 

(ii) If the line through q' and (0, ^min) does not intersect the interior of U{p), then we return 
p with elevation Zuim, unless Zmm > high{p). 

(iii) Otherwise, we conduct a binary search on Z{p) as indicated above to find the lowest 
intersection (0, z) of the vertical axis and a tangent of U{p) through q'. If z > high{p), we 
do not include p in the result, otherwise, we return p with elevation max(z, low{p)). Note 
that we do not need to remove itself from U{p) (and Z{p)) in this procedure, since it 
will never lose if it competes with itself. 

The computations can be done in time logarithmic in the degree of p. □ 

4.1.4 Correctness and running time of the complete algorithm 

Theorem 2 After precomputations in 0(n log n) time and 0{n) space, the algorithm 
PotentialWS((5) computes the potential watershed Wy {Q) of a set of nodes Q and its canon- 
ical realization RuiQ) in time 0(n log n), where n is the number of edges in the terrain. 

Proof: The algorithm searches the graph starting from the nodes of Q. At each point in time 
we have three types of nodes. Nodes that have been extracted from the priority queue have a 
finalized elevation, a node that is currently in the priority queue but was never extracted (yet) 
has a tentative elevation, other nodes have not been reached. 

We show that when (p, z) is first extracted from the priority queue in Algorithm [T| p is 
indeed contained in the potential watershed of Q, and the elevation z is the lowest possible 
elevation of p such that water flows from p to some node in Q. To this end we use an induction 
on the points extracted, in the order in which they are extracted for the first time. 

The induction hypothesis consists of two parts: 

(i) There exists a realization R and q & Q such that elevpiip) = and R induces a fiow path 
TT from p to q which only visits vertices that have been extracted from the priority queue. 

(ii) There exists no realization R and q & Q such that elevji{p) < z and Pj^q- 

If a node q G Q is extracted with z = low{q), then the claims hold trivially. Note that the 
first extraction from the priority queue must be of this type. 

If p is extracted from the priority queue for the first time and p ^ Q, then there must be 
at least one node p' that was extracted earlier, such that Expand (p', z'), for some elevation z' , 
resulted in p having the tentative elevation z. By induction, there exists a realization R' and 
q £ Q, such that elevR'ip') = z' ., there is a fiow path tt from p' to q in i?', and vr does not 
include p. 
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To see (i), we construct a realization R by modifying R' as follows: we set elevji{p) = z, 
and we set elevfi{r) = high{r) for each neighbor r of p that does not lie on vr. In comparison to 
R' , only p and its neighbors may have a different elevation in R. Since elevR^p) = z > z' is still 
at least as high as the elevation of any node on vr, water will still flow along the path vr from 
p' to q. By the definition of Expand, none of the neighbors of p that are set at their highest 
elevation can out-compete p' as a steepest-descent neighbour of p. Therefore, the steepest- 
descent neighbour of p in R' must be one of the nodes on vr. Thus, water from p will flow onto 
vr, and thus, to q. 

Next we show (ii). Suppose, for the sake of contradiction, 
there is a realization R such that elevji{p) < z and there is a 
flow path from p to a node q G Q. Consider two consecutive 
nodes r and s on this path, such that r has not been extracted 
before but s has been previously extracted (it may be that r = p 
and/or s G Q). Note that flow paths have to be monotone in the 
elevation. We argue that this path cannot stay below z in any 
realization. Since r is a neighbor of s, it has been added to the 

priority queue during the expansion of s. Let the tentative elevation of r that resulted from this 
expansion be z^. By induction, since the elevation of s is finalized, Zr is a lower bound on the 
elevation of r for any flow path that follows the edge (r, s) and then continues to a node in Q 
in any realization. However, Zr > z, since r was not extracted from the priority queue before p. 
Therefore, a path from p to q that contains r with elevji{p) < z cannot exist. This proves (ii). 

It follows that the algorithm outputs all nodes of (Q) together with their elevations in 
RuiQ). 

As for the running time, computing and storing U {p) and Z{p) for a node p of degree d takes 
0(dlog(i) time and 0{d) space. Since the sum of all node degrees is 2n, computing and storing 
U{p) and Z{p) for all nodes p thus takes 0(n log dmax) time and 0{n) space in total, where dmax 
is the maximum node degree in the terrain. While running algorithm PotentialWS(Q), each 
node is expanded at most once. By Lemma [5| Expand (g', z') on a node q' of degree d takes 
time 0(dlog dmax)- Thus, again using that all nodes together have total degree 2n, the total 
time spent on expanding is 0(n log dmax) = 0{n\ogn). Each extraction from the priority queue 
takes time O(logn) and there are at most 0{n) nodes to extract. Therefore PotentialWS 
takes time 0(n log n) overall. □ 

For grid terrains, dmax = 0(1), and thus, the slope diagram computations take only 0(1) 
time per expansion. In fact, since we only need to expand nodes that are in Wy (Q), we could 
actually compute Wy (Q) in 0{k\ogk) time, where k = IWy (Q) |. Alternatively, we can use 
the techniques from Henzinger et al. [15] for shortest paths to overcome the priority queue 
bottleneck, and obtain the following result (details in Appendix [A]) : 

Theorem 3 The canonical realization of the potential watershed of a set of cells Q in an im- 
precise grid terrain of n cells can be computed in 0{n) time. 



4.2 Potential downstream areas 

Similar to the potential watershed of a set Q, we can define the set of points that potentially 
receive water from a node in Q. Let 
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Naturally, a canonical realization for this set does not necessarily exist, however, it can be 



computed in a similar way as described in Section 4.1 using a priority queue that processes 
nodes in decreasing order of their maximal elevation such that they could still receive water 
from a node in Q. The algorithm is the same as Algorithm [T| except that in the first line the 
nodes are enqueued with their highest possible elevation, in line 3 we dequeue the current node 
with the largest key and we use the following subroutine in line 6. 

Definition 3 Let ExPANDDowN(g', z') denote a function that returns for a node q' and an 
elevation z' G [low{q'), high{q')] a set of pairs of nodes and elevations, which includes the pair 
{p,z) if and only if p £ N{q'), there is a realization R with elevii[q') € [low{(^), z'] such that 
q' ^p, and z is the maximum elevation of p over all such realizations R. 




Lemma 6 We can compute the function ExPANDDowN(g', z') in 0{d\ogd) time, where d is 
the node degree of q' . 



Proof: Consider the slope diagram of q' as defined in Section 4.1.3 Let zq be mmhigh{p) 
over all neighbours p of q'; note that this is the vertical coordinate of the lowermost point of 
U{q'). Let q' = (0, z') and consider its lower tangent to U{q'). Let pi be the corner of U{q') 
that intersects the tangent. Similarly, let pj be the corner of U{(^) that intersects the tangent 
through (0, max(/oio(g'), zo)). Let W{q') be the intersection of the halfplanes above these two 



tangents and the halfplanes Hi, . . . ,Hj as defined in Section 4.1.3 Clearly, a neighbor of q' can 
have a steepest descent edge from q' , for some elevation of q' in [low{p), z'], if and only if its 
representative in the slope diagram, lies below W{q') or on the boundary of W{q'). To compute 
the neighbors of q' and their elevations as they should be returned by ExPANDDowN(g', z'), 
we test each neighbor p of q' as follows. We find the point p' = {\pq'\,z) that is the projection 
from p down onto the boundary of W{q'). If z > low{p), we return {p, z), otherwise we do not 
include p in the result. 

The slope diagram with W{q') can be computed 0{dlogd) time. The neighbors p of q' can 
be sorted by increasing distance from q' in the xy-projection in 0{dlogd) time; after that, the 
projections of all points p can be computed in 0{d) time in total by handling them in order of 
increasing distance from q' and walking along the boundary of W{q') simultaneously. □ 

Theorem 4 Given a set of nodes Q of an imprecise terrain, we can compute the set Du(Q) 
time 0(n log n), where n is the number of edges in the terrain. 
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Proof: The algorithm searches the graph starting from the nodes of Q. As in the algorithm for 
potential watersheds, nodes that have been extracted from the priority queue have a finalized 
elevation; nodes that are currently in the priority queue but were never extracted (yet) have 
tentative elevations. However, this time these elevations are not to be understood as elevations 
of the nodes in a single realization, but simply as the highest known elevations so that the nodes 
may be reached from Q. 

The induction hypothesis is symmetric to the hypothesis used for potential watersheds: we 
show that when {p, z) is first extracted from the priority queue, p is indeed contained in the 
potential downstream area of Q, and the elevation z is the highest possible elevation of p such 
that water flows from some node in Q to p. Again, the induction is on the points extracted, in 
the order in which they are extracted for the first time. 

The induction hypothesis consists of two parts: 

(i) There exists a realization R and q & Q such that elevji{p) = z, there is a flow path vr from 
g to p in R, and tt only visits vertices that have been extracted from the priority queue. 

(ii) There exists no realization R and q G Q such that elevR{p) > z and q^P- 

If a node g G Q is extracted with z = high{q), then the claims hold trivially. Note that the 
first extraction from the priority queue must be of this type. 

If p is extracted from the priority queue for the first time and p ^ Q, then there must be at 
least one node p' that was extracted earlier, such that ExPANDDowN(y, z'), for some elevation 
z' , resulted in p having the tentative elevation z. By induction, there exists a realization R' 
and q £ Q, such that elev^i{p') = z' , there is a flow path vr from q to p' in R' , and vr does not 
include p. 

So far the proof is basically symmetric to that of Theorem [2] However, to see (i), we 
need a different construction. Let z" < z' be the elevation such that water flows from p' to 
p in the realization R" with elevjiii{p') = z" , elevjiii{p) = z, and elevfi"{p") = high{p") for 
all other nodes p" . Note that z" exists by definition of ExpandDown. We now construct a 
realization R by modifying R' as follows: we set elevR{p') = z", we set elevR{p) = z, and we 
set elevii{r) = high(r) for each neighbor r of p' such that r ^ p and r does not lie on tt. In 
comparison to R' , only two nodes in R may have lower elevation, namely p and p' . Therefore, 
water will still flow along the path vr from q until it either reaches p' , or a vertex that now has 
p or p' as a new steepest-descent neighbour. Thus, in any case, there is a flow path from q to 
either p or p' . If the flow path reaches p' , then, by definition of ExpandDown, none of the 
neighbours of p' that are set at their highest elevation can out-compete p as a steepest-descent 
neighbour oi p' . Of course, the neighbours of p' that lie on vr cannot out-compete p either, since 
these neighbours have elevation at least as high as p' . Therefore, p must be a steepest-descent 
neighbour of p' in R' , and water from p' will flow to p. Thus, in any case, water from q will 
reach p in R' along a path that is a prefix of vr, followed by an edge to p. This proves part (i) 
of the induction hypothesis. 

The proof of part (ii) is completely analogous to the proof of Theorem [2| 
It follows that the algorithm outputs all nodes of 'D\j(Q). The running time analysis is 
analogous to Theorem [2| □ 

4.3 Persistent watersheds 

In this section we will give a definition of a minimal watersheds, and explain how to compute 
it. Recall that the potential (maximal) watershed of a node set Q is defined as the set of nodes 
that have some flow path to a node in Q. We can write this as follows 

Wu (Q) = {p : 3 vr e U{31t),tt 3 p 3 q £ Q : p n^q} . 
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"potential local minimum" 
1 

Figure 9: An example of a 1.5 dimensional imprecise terrain, where the core watershed can be 
arbitrarily reduced by oversampling. The node p cannot be in the core watershed of any other 
node. 



An analogous definition to this would be 

Wn((3) = {p: V7rGn(3iT),vr9p3gGQ:p^g}. 

This is the set of nodes p from which water flows to Q via any induced flow path that contains p. 
We call this the core watershed of Q. 

However, this definition seems a bit too restrictive. Consider the case of a measuring device 
with a constant elevation error, used to sample points in a gently descending valley. It is 
possible that, by increasing the density of measurement points, we can create a region in which 
imprecision intervals of neighbouring nodes overlap in the vertical dimension, and thus each 
node could become a local minimum in some realization. Thus, water flowing down the valley 
could, theoretically, "get stuck" at any point, and thus, the minimum watershed of the point q 
at the bottom of the valley would contain nothing but q itself, see Figure |9j 

Nevertheless, it seems clear that any water flowing in the valley must eventually reach q 
(possibly after flooding some local minima in the valley), since the water has nowhere else to 
go. This leads to an alternative definition of a minimal watershed, after we rewrite the definition 
of the core watershed slightly. Observe that the following holds for the complement of the core 
watershed. 

(Wn(Q))^ = {p: 3 TT G n(3iT),vr 9 p -3 g G Q : p^^^} 

Thus, the core watershed of Q is the complement of the set of nodes p, for which it is possible 
that water follows a flow path from p that does not lead to Q. Assume there exists a suitable 
set of alternative destinations S, such that we can rewrite the above equation as follows 

(Wn(Q))' ={p: 3 ^ G U{JIt), tt 3 p 3 s e S : {p jr^s) A {tt[p, s] n Q = 0)} . 

Note that the right hand side is equivalent to the set 

V^}?{S):= U [J{p:ipji,s)Ai7T[p,s]nQ = (l>)} (1) 
7ren(3?T) se5 

We call the set in Equation [T] the Q- avoiding potential watershed of a set of nodes S and 
we denote it with w}f^ (S). This is the set of nodes that have a potential flow path to a node 

^Interestingly, there are some parallels to observations made in the GIS literature. Firstly, Hebeler et al. [ij 
observe that the watershed is more sensitive to elevation error in "flatlands" . Secondly, simulations have shown 
that also potential local minima or "small sub-basins" can severely affect the outcome of hydrological computa- 
tions [18j . 
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s S 5 that does not pass through a node of Q before reaching s. Note that it is still possible for 
those flow paths to intersect Q, as long as this happens outside the subpath between p and s. 

It remains to identify the set of alternative destinations S. Since every flow path ends in 
a local minimum, the set of potential local minima clearly serves as such a set of destinations. 
Let V^-^ be the union of all sets r such that there exists a realization in which all nodes of r 
have the same elevation, r is a local minimum, and rflQ = 0. However, it is also safe to include 
the nodes that do not have any flow path to Q, which is the complement of the set Wy (Q). It 
follows for the core watershed: 

Wn(Q) = ( ( U ( Wu {Q)r ) y 
Note that we can rewrite this as follows: 

Wn(Q) = ( W^'^ ( ( Wu {Q)f ) \ W^'^ ( lAO n Wu (Q) ) 

Based on the above considerations we suggest the following alternative definition of a mini- 
mal watershed. 

Definition 4 The persistent watershed of a set of nodes Q is defined as 

Wn(Q):=(w\«( (Wu (Q)r ) 



The shaded area in Figure [9] indicates what would be the persistent watershed of q in this case: 
these are the nodes that can never be high enough so that water from those nodes could escape 
from the potential watershed of q. 

To compute the persistent watershed efficiently, all we need are efficient algorithms to com- 
pute potential watersheds and Q-avoiding potential watersheds. We have already seen how to 



compute Wu (Q) efficiently in Section 4.1 Note that the Q-avoiding potential watershed of S 
is different from the potential watershed of S in the terrain T' that is obtained by removing 
the nodes Q and their incident edges from T. The next lemma states that we can also compute 
Q-avoiding potential watersheds efficiently. 

Lemma 7 There is an algorithm which outputs the Q-avoiding potential watershed of S and 
takes time 0(n log n), where n is the number of edges of the terrain. 

Proof: We modify the algorithm to compute the potential watershed of S as shown in Algo- 
rithm [l| such that, each time the algorithm extracts a node from the priority queue, this node 
is discarded if it is contained in Q. Instead, the algorithm continues with the next node from 
the priority queue. Clearly, this algorithm does not follow any potential flow paths that flow 
through Q. However, the nodes of Q are still being considered by the neighbors of its neighbors 
as a node they have to compete against for being the steepest descent neighbor. It is easy 
to verify that the proof of Theorem [2] also holds for the computation of Q-avoiding potential 
watersheds. □ 

By applying Theorem [2] and Lemma [7j we obtain: 

Theorem 5 We can compute the persistent watershed Wfri(Q) of Q in time 0(n log n), where 
n is the number of edges of the terrain. 
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5 Regular terrains 



We extend the results on imprecise watersheds in the network model for a certain class of 
imprecise terrains, which we call "regular". We will first define this class and characterize it. 
To this end we will introduce the notion of imprecise minima (see Definition [s]) , which are the 
"stable" minima of an imprecise terrain, regular or non-regular. In Section 5.2 we will describe 
how to compute these minima and how to turn a non-regular terrain into a regular terrain. 
In the remaining sections, we discuss nesting properties and fuzzy boundaries of imprecise 
watersheds. Furthermore, we observe that regular terrains have a well-behaved ridge structure, 
that delineates the main watersheds. 

The main focus of this section is on the extension of the results in section Section |4l Some 
of the concepts introduced here could also be applied to the surface model, however, we confine 
our discussion to the network model. 



5.1 Characterization of regular terrains 

We first give a definition of a proper minimum in an imprecise terrain. 



Definition 5 A set of nodes S in an imprecise terrain T is an imprecise minimum if S 

contains a local minimum in every realization of T, and no proper subset of S has this property. 



Now a regular imprecise terrain is defined as follows: 



Definition 6 An imprecise terrain T is a regular imprecise terrain, if every local minimum 
of the lowermost realization R~ ofT is an imprecise minimum ofT. 

Any imprecise minimum S* on a regular terrain is a minimum in R~ . Indeed, if S would 
not be a minimum on , then, by Definition [sj it would contain a proper subset S' that is 
a minimum on R~ while S' is not an imprecise minimum of T - but this would contradict 
Definition [6} Now, we observe: 

Observation 1 Let S be an imprecise minimum on a regular terrain. Then each node s £ S 
has the same elevation lower bound low(s). Furthermore, for each subset S' <Z S we have 
Wu (5') = Wu (5) and Wr- {S') = Wr- (S). 

We derive a characterization of imprecise minima. For this, we introduce proxies. 

Definition 7 A proxy of an imprecise minimum S is a node p £ S, such that there are no 
realizations R and nodes q ^ S such that Pj^q. 

Thus, water that arrives in a proxy of an imprecise minimum S, can never leave S anymore. 
This implies that the proxy is not in the potential watershed of any set of nodes that lies entirely 
outside S. The following lemma states that every imprecise minimum contains a proxy. 

Lemma 8 Let the bar of a set S be bar{S) = minxes' high{s). A set S is an imprecise minimum 
if and only if (i) bar[S) < min^g^^g) low{t) and (ii) no proper subset S' of S has this property. 
Every imprecise minimum has a proxy. 

Proof: First observe that condition (i) implies that S contains a local minimum in any realiza- 
tion. 
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If S is an imprecise minimum, then, by definition, it contains a local minimum in any 
realization and no proper subset S' of S has this property. We argue that this implies (i) and 
(ii) for S. 

To prove (i), consider the following realization R: For all nodes r E 5" we set elevji{r) = 
max{bar{S), low{r)), and for all nodes t G A^(-S') we set elevji{t) = low{t). Now suppose, for 
the sake of contradiction, that there is a proper subset S' of S that is a local minimum in 
R. Like all nodes of S, the local minimum S' must have elevation at least bar{S); each node 
t G N{S') must be set at a higher elevation low{t). If we would remove the nodes of N{S') from 
5, the imprecise minimum S would be separated into several components, including at least one 
component S" that contains a node s with high{s) = har{S). This component S" is a proper 
subset of S. Its neighbourhood N{S") consists of nodes from N[S) and N{S'), all of which 
have an elevation lower bound strictly above bar{S) = min^g^" high{s). Thus S" C S meets 
condition (i) and contains a local minimum in any realization, contradicting the assumption 
that S is an imprecise minimum. If follows that no proper subset S' of S is a local minimum in 
R; therefore S must be a local minimum as a whole, which implies (i). 

To prove (ii), assume, for the sake of contradiction, that S contains a proper subset S' such 
that har{S') < min(gjv(S") low{t). Thus, S' would contain a local minimum in any realization, 
and S would not be an imprecise minimum; hence (ii) must hold for S. 

Now we argue that, if (i) and (ii) are met, then S is an imprecise minimum. Recall that if 
condition (i) is met, then S contains a local minimum in any realization. Now assume, for the 
sake of contradiction, that there exists a proper subset S' that always contains a local minimum. 
Let S' be a smallest such subset of S. We have that S' is an imprecise minimum, and therefore, 
as we proved above, it holds that bar{S') < mm^^]\f(^gi-j low{t), which contradicts that condition 
(ii) holds for S. Hence, there is no proper subset S' of S that always contains a local minimum; 
therefore S is an imprecise minimum. 

As a proxy of an imprecise minimum S, we take any node s such that high(s) < min^gjv(5) low{t). 
By part (i) of the lemma, such a node s always exists. Since s lies below any node of N{S) in 
any realization, there are no realizations R and nodes q ^ S such that Sj^q; thus s is a proxy 
of 5. □ 

5.2 Computing proxies and regular terrains 

Any imprecise terrain can be turned into a regular imprecise terrain by raising the lower bounds 
on the elevations such that local minima that violate the regularity condition are removed 
from R~ . Indeed, in hydrological applications it is common practice to preprocess terrains by 
removing local minima before doing flow computations [25j. To do so while still respecting 
the given upper bounds on the elevations, we can make use of the algorithm from Gray et 
al. [11] . The original goal of this algorithm is to compute a realization of a surface model that 
minimizes the number of local minima in the realization, but the algorithm can also be applied 
to a network model. It can easily be modified to output a proxy for each imprecise minimum of 
a terrain. Moreover, the realization M computed by the algorithm has the following convenient 
property: if we change the imprecise terrain by setting low{v) to elevM{v) for each node, we 
obtain a regular imprecise terrain. 

The algorithm The algorithm proceeds as follows. We will sweep a horizontal plane upwards. 
During the sweep, any node is in one of three states. Initially, each node is undiscovered. Once 
the sweep plane reaches low{v), the state of the node changes to pending. Pending nodes are 
considered to be at the level of the sweep plane, but they may still be raised further. During the 
sweep, we will always maintain the connected components of the graph induced by the nodes 
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that are currently pending; we call this graph Gp. As soon as it becomes clear that a node 
cannot be raised further or does not need to be raised further, its final elevation on or below 
the sweep plane is decided and the node becomes final. More precisely, the algorithm is driven 
by two types of events: we may reach low{v) for some node v, or we may reach high{v) for some 
node V. These events are handled in order of increasing elevation; Zow;(t')-events are handled 
before high{v)-events at the same elevation. The events are handled as follows: 

• reaching low^v): we make v pending, and find the component S of Gp that contains v. If 
V has a neighbour that is final, we make all nodes of S final at elevation low{v). 

• reaching high{v): if v is final, nothing happens; otherwise we report v as a. proxy, we find 
the connected component S of Gp that contains v, and we make all nodes of S final at 



Gray et al. explain how to implement the algorithm to run in 0(n log n) time, jllj. 

Lemma 9 Given an imprecise terrain T, (i) all nodes reported by the above algorithm are 
proxies of imprecise minima, and (ii) the algorithm reports exactly one proxy of each imprecise 
minimum ofT. 

Proof: We first prove the second part, and then the first part of the lemma. 

(ii) Let S be an imprecise minimum. Let v be the node in S which was the first to have its 
high{v)-event processed. By Lemmapl u is a proxy of S and we have high{v) < ?7T.mjg7v(s) low{t). 
Hence, when high{v) is processed, the component of Gp that contains v does not contain any 
nodes outside S, and the high{v)-event is the first event to make any nodes in this component 
final. Thus, v is reported as a proxy. Furthermore, no node s £ S can have low{s) > high{v), 
otherwise bar(S \ {s}) = high{v) < min^fz^s, N{s)} low{t) < min^fzN (^s\{s}) low{t), and thus, by 
Lemma [sj S would not be an imprecise minimum. Hence, when the high{v)-event is about to 
be processed, all nodes of S have been discovered and are currently pending. The high{v)-event 
makes all nodes of S final; thus, any high{s)-events for other nodes s £ S will remain without 
effect and no more proxies of S will be reported. 

(i) Let u be a node that is reported as a proxy in a high{v)-event. We claim that the con- 
nected component 5 of Gp that contains v at that time, is an imprecise minimum. Indeed, 
by definition of Gp, all nodes of S are pending, and thus high{v) = mms,=s high{s) = bar{S). 
Furthermore, because 5 is a connected component of G(P), all nodes t £ N{S) must be ei- 
ther undiscovered or final. In fact, the algorithm maintains the invariant that no neighbor of a 
finalized node is pending; since all nodes in S are pending, all nodes t G N{S) must be undiscov- 
ered. Therefore high{v) < min^gjv(5) low{t). Because all low {t)-events at the same elevation as 
high{v) are processed before the high{v)-event is processed, we actually have a strict inequality: 
high{v) < minig7v(s) low{t). It follows that S satisfies condition (i) of Lemma [sj Furthermore, 
no proper subset S' of S has this property, otherwise, by the analysis given above, a proxy for 
S' would have been reported already and the nodes from S' would have been removed from Gp 
at that time. Hence, S also satisfies condition (ii) of Lemmajsj and S is an imprecise minimum, 
with V as a proxy. □ 

Lemma 10 Let M be the realization of a terrain T as computed by the algorithm described 
above. Let T' be the imprecise terrain that is obtained from T by setting low{v) = elevM{v) for 
each vertex v. The terrain T' is a regular imprecise terrain. 

^This is a small variation: the algorithm as described originally by Gray et al. would make the elevations 
final at highiv) = minsgs high{s). However, in the current context we prefer to make the elevations final at 
maxsgs low{s), to maintain as much of the imprecision in the original imprecise terrain as possible. 
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Proof: Note that M is the lowermost reaUzation of T' . Consider any local minimum S of M. 
Observe that the algorithm cannot have finalized the elevations of the last pending vertices of S 
in a low{v)-event, because then we would have v G S and v must have a neighbor t ^ S that was 
finalized before v; hence elevuit) < elevM{v) and S would not be a local minimum. Therefore, 
the algorithm must have finalized the last elevations of the vertices of /S in a high{v)-event 
for a vertex v £ S. Furthermore, each vertex t G N{S) must have been undiscovered at that 
time; otherwise t would have become part of the same component as the vertices of S before 
its elevations were finalized, or t would have been finalized before v: in both cases S would not 
be a local minimum. Hence we have low{t) > high{v) for each vertex t G N{S), and thus, S 
is a local minimum in every realization of T or T' . Furthermore, no proper subset of S' of S 
contains a local minimum in every realization of T' , since in particular, in M the set S is a local 
minimum and therefore no proper subset S' of 5 is a local minimum. Thus, by Definition [5] and 
Definition [6| T' is a regular terrain. □ 

5.3 Nesting properties of imprecise watersheds 

To be able to design data structures that store imprecise watersheds and answer queries about 
the flow of water between nodes efficiently, it would be convenient if the watersheds satisfy the 
following nesting condition: if p is contained in the watershed of q, then the watershed of 
p is contained in the watershed of q. Clearly, potential watersheds do not satisfy this nesting 
condition, while core watersheds do. However, in general, persistent watersheds, are not nested 
in this way. We give a counter-example that uses a non-regular terrain in the next lemma before 
proving the nesting condition for persistent watersheds in regular terrains later in this section. 



Lemma 11 There exists an imprecise terrain with two nodes p and 1 



such that p G y^(^{q) and 



Proof: We give an example of a non-regular terrain that has this property. Refer to Figure 10 
The persistent watershed of p as shown in red is not completely contained in the persistent 
watershed of q as shown in blue. The left figure gives a top-view. All edges have unit length, 
except for the edge between w and q. The right figure shows the fixed elevations of s,t,t',u,v 
and w, the elevation intervals of p, q and r, and the correct horizontal distances on all edges 
except \pv\ and \qv\. The red outline delimits Wy (p) = {p,q,r, s,t,v,w}. The red dashed 
outline delimits Wfr|(p) = {p, s, v}. The blue outline delimits Wy (q) = {p, q, s, v, w}. The blue 
dashed outline delimits Wfr|(g) = {p,q,v}. □ 
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Figure 10: Counterexample of a non-regular terrain to the nesting condition of persistent wa- 
tersheds. 
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The following lemmas will prove that on regular imprecise terrains persistent watersheds do 
satisfy the nesting condition. 



Lemma 12 Let Q be a set of nodes in a regular imprecise terrain, then Wfr|(Q) C Wj^-{Q). 

Proof: Consider a flow path from a node p G y^i^{Q) in . By the definition of persistent 
watersheds, the path cannot leave Wy (Q) without going through Q. 

If the path reaches a local minimum S without going through any node of Q, we claim that 
this local minimum must contain a node q £ Q. To prove this claim, assume, for the sake of 
contradiction, that there is a Q-avoiding flow path to a local minimum S such that Q Ci S = f/i. 
Thus there would be a Q-avoiding flow path to any node of S, and in particular, to a proxy 
s € 5; because the terrain is regular, S must be an imprecise minimum (by Definition [g]) , and 
therefore a proxy is guaranteed to exist by Lemma[8| As observed above, s is not in the potential 
watershed of any set of nodes outside S; in particular, s is not in Wy (Q). This implies that 
there is a realization in which a flow path from p leaves Wy (Q) without going through any 
node of Q, contradicting the assumption that p £ W(t^{Q). Hence, if a flow path from p in R~ 
reaches a local minimum S without going through any node of Q, then S must contain a node 
q £ Q, and there would also be a flow path from p to q. 

Therefore, from any node p € "Wfr|(Q) there must be a flow path to a node q £ Q in R~ , and 
thus, Wn(Q) C W^-(Q). □ 

Lemma 13 Let Q be a set of nodes on an imprecise terrain, and let P C W^-(Q). Then 
Wu (P) C Wu (Q). 

Proof: Let R be the watershed-overlay of ^^r^ (-P) (-^) ^R- (Q) • Consider a node r £ Wy (P) 
and a flow path vr from r to a node p £ P in Ru{P). Let vr' be the maximal prefix of vr such 
that the nodes of vr' have the same elevation in i?u(P) and R, and let vr" be the maximal prefix 
of tt' such that vr" is still a fiow path in R. We distinguish three cases: 

• If vr' = vr" is empty, then r has lower elevation in Ru{P) than in R, so r must be in 
Wr-(Q). 

• If vr' = vr" = TT, then fiow from r reaches a node p £ P CI W^- (Q) in R. 

• Otherwise, let {u,v) be the edge of tt such that u is the last node of vr". Now v is not 
on vr", so in R, flow from u either still follows {u,v) but elev^{v) < elev r^(^p~^{v) , or flow 
from u is diverted over an edge (u, to another node v with elev^iv) < elev ji^(^p){v) . In 
either case, from u we follow an edge to a node of which the elevation in R is lower than 
in Ru{P); therefore this must be a node of "W^-(Q). 

In any case, there is a flow path from r to a node of Wj:j-(Q). From here, there must a path to 
a node q £ Q, since every flow path within W^- (Q) in R~ is also a flow path in R. Thus there 
is flow path from r to g in R, and thus, r £ Wy (Q). This proves the lemma. □ 

Lemma 14 (persistent watersheds are nested) Let Q be a set of nodes in a regular imprecise 
terrain, and let P C W(^{Q). Then W^{P) C Wr(Q). 

Proof: Assume for the sake of contradiction that there exists a node s £ 'Wfri(P), such that 



s ^ Wfr|((5). Clearly, s £ Wy (P) and by Lemma 13 and Lemma 12, it holds that s £ Wy (Q). 



Furthermore, s must have a flow path vr to a point r ^ Wy (Q), which does not pass through 



a node of Q, refer to Figure 11 By Lemma 13 and Lemma 12, Wy (P) C Wy (Q), and thus 



r ^ Wu (P). Furthermore, since s £ W(r^(P), the subpath 7r[s,r] must include a node p £ P. 
This contradicts the fact that p £ W^{Q), since vr[p, r] n Q = 0. □ 
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Figure 11: Persistent watersheds are nested on regular terrains. Illustration to the proof of 
Lemma [TH 



5.4 Fuzzy watershed boundaries 



Lemma [12] and Lemma 13 also allow us to compute the difference between the potential and a 
persistent watershed of a set of nodes Q efficiently, given only the boundary of the watershed 
of Q on the lowermost realization of the terrain. We first define these concepts more precisely. 

Definition 8 Given a realization R, and a set of nodes Q, let Xr{Q) he the directed set of 
edges {u,v) such that u G W/j((5) and v ^ W/j(Q). We call Xr{Q) the watershed boundary 
of Q in R. Likewise, we define the fuzzy watershed boundary of Q as the directed set of 
edges {u,v) such that u € Wy (Q) and v ^ Wfr|((5) and we denote it with X[j(Q). We call the 
set Wu (Q) \ "WfT|((5) the uncertainty area of this boundary. 

We will now discuss how we can compute the uncertainty area of any fuzzy watershed 
boundary efficiently. 

Algorithm to compute the uncertainty area of a watershed. Assume we are given 
Xr-{Q)- We will compute the set Wy (Q) \ 'Wi^iQ) with Algorithm [T| modified as follows. 
Instead of initializing the priority queue with the nodes of Q, we initialize in the following 
way. For each edge {u, v) G Xji- (Q), we use the slope diagram of u to determine the minimum 
elevation Zu of u, such that there is a realization in which water flows on the edge from u to 
V. If there exists such an elevation Zu, we enqueue u with elevation (and key) Zu- Similarly, we 
use the slope diagram of v to determine the minimum elevation z^ of v, such that water may 
flow on the edge from v to u. If z^ exists, we enqueue v with elevation (and key) z^. After 
initializing the priority queue in this way, we run Algorithm [T] as written. 

Lemma 15 // the terrain is regular, the algorithm described above computes Wy (Q) \ 'VVfTi(Q). 
This is the uncertainty area of the fuzzy watershed boundary of Q. 

Proof: Observe, following the proof of Theorem [2j that for any node p output by the above 
algorithm, there are a realization Rg and a node s which was in the initial queue with elevation 
Zs, such that elevR^{s) = Zs and Rs induces a flow path vr from p to s. Let {u,v) be the edge 
of X^-[Q) which led to the insertion of s G {u,v} into Q with elevation Zg- Let t be the other 
node of {u,v), that is, t = {u,v} \ {s}. let Rt be the realization obtained from Rs by setting 
elevR^{t) = low{t). Observe that, by our choice of Zs, the realization Rt now induces a flow 
path from p to t. We will now argue that (i) p G Wy (Q), and (ii) p ^ 'WfT|(Q). 

(i) The existence of Ru implies that p G Wy (u); since u G W^- {Q) (by deflnition of X^- {Q) 



this implies p G Wy (Q) (by Lemma 13). 



(ii) By deflnition of Xji-{Q), there is no flow path from v to Q on R . Hence, any flow 
path from v on R~ must lead to a local minimum S that does not contain any node of Q, and 



26 



by Definition |6j each such local minimum S is an imprecise minimum. Now, by Lemma [8j each 
such local minimum S contains a proxy s, which is, by Definition [7| not contained in Wy {Q). 
Thus there is a flow path from v that does not go through any nodes of Q and leads to a proxy 
s ^ Wu (Q). Hence, by Definition [4| p ^ W^{Q). 

Next, we will argue that if p £ Wy (Q) and p ^ "WfT|(Q), the algorithm will output p. We 
distinguish two cases. 

If p G W_R-(Q); then, because p ^ y^piiQ), there must be a flow path on from p to a 
minimum S that does not contain any node of Q. By Definition [6j Lemma [8] and Definition [7| 
there will then be a flow path from p to a proxy s £ S that lies outside Wy (Q), and thus, 
outside W^- (Q). 

If p ^ y^R-iQ), then, because p £ "Wy (Q), there must be a realization in which there is a 
flow path from p to Q, and thus, from p to W^- (Q). 

In both cases, there is a realization in which there is a flow path from p that traverses an 
edge (u, v) E Xj^- (Q), either from u to v or from v to u. The algorithm reports at least all such 
points p. 

This completes the proof of the lemma. □ 

Note that if all nodes have degree 0(1), the running time of the above algorithm is linear 
in the size of the input (Af/j- (Q)) and the output (Wu (Q) \ 'Wfri(Q)). When a data structure is 
given that stores the boundaries of watersheds on so that they can be retrieved efficiently, 
and the imprecision is not too high, this would enable us to compute the boundaries and sizes of 
potential and persistent watersheds much faster than by computing them (or their complements) 
node by node with Algorithm [TJ 

We can use the same idea as above to compute an uncertain area of the watershed boundaries 
between a set of nodes Q. More precisely, given a collection of nodes Q such that no node q £ Q 
is contained in the potential watershed of another node q' £ Q, we can compute the nodes that 
are in the potential watersheds of multiple nodes from Q. 

Algorithm to compute the uncertainty area between watersheds. Let Q be {qi, qk} 
and let G' be the graph induced by the potential watershed of Q. The algorithm is essentially 
the same as algorithm that computes the uncertainty area of a single watershed's boundary — 
the main difference is that now we have to start it with a suitable set of edges X on the fuzzy 
boundaries between the watersheds of the nodes of Q. More precisely, X should be an edge 
separator set of G' , which separates the nodes of G' into k components G'i,...,G'f^ such that 
nodes of each component G'^ are completely contained in Wy (qi). 

We obtain X with the following modification of Algorithm [TJ For each node p we will 
maintain, in addition to a tentative elevation z, a tentative tag that identifies a node q £ Q 
such that there is a realization R with elevR{p) = z and P^q- We initialize the priority 
queue of Algorithm [l] with all nodes q £ Q, each with tentative elevation low{z) and each 
tagged with itself. The first time any particular node q' is extracted from the priority queue, 
we obtain not only its final elevation but also its final tag q from the queue, and each pair 
{p, z) £ Expand((7', z') is enqueued with that same tag q. At the end of Algorithm[lJ we obtain 
the set of nodes in Wy {Q) together with their elevations in the canonical realization R\j{Q) 
and with tags, such that any set of nodes tagged with the same tag q £ Q forms a connected 
subset of (g). We now extract the separator set X by identifying the edges between nodes 
of different tags. 

Having obtained X, we compute the union of the pairwise intersections of the potential 
watersheds of qi,...,qk as follows. Again, we use Algorithm [TJ This time the priority queue is 
initialized as follows. For each edge {u, v) £ X, we use the slope diagram of u to determine the 
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minimum elevation of u, such that there is a reahzation R with elevR{v) = elev ji^(^q-j{v) in 
which water flows on the edge from u to v. If there exists such an elevation Zu, we enqueue u 
with elevation (and key) Zu- Similarly, we use the slope diagram of v to determine the minimum 
elevation z^ of v, such that water may flow on the edge from v to u at elevation elev r^(^q){u) . 
If Zy exists, we enqueue v with elevation (and key) z^. After initializing the priority queue in 
this way, we run Algorithm [T] as written, and output the result. 

Lemma 16 Given a set of nodes qi, . . . ,qk of an imprecise terrain, such that Qi ^ Wy (qj) for 
any i ^ j and 1 < i,j < k, we can compute the set UiUjyjC^u {Qi) H (qj)) in O(nlogn) 
time, where n is the number of edges of the imprecise terrain. 

Proof: The separator set X is obtained in 0(n log n) time by running the modified version of 
Algorithm [T] and one scan over the graph to identify edges between nodes with different tags. 
Computing the union of the pairwise intersections of the potential watersheds of qi, ...,qk with 
the modified Algorithm [l] takes O(nlogn) time again. 

By the same arguments as in the proof of Lemma [Tsj we can observe the following: for 
any node p output by the above algorithm, there is an edge {u, v) £ X, a realization Ru with 
elevji^{u) = elev ji^(^Q-){u) and Pj^u, and a realization R^ with elevji^{v) = eZeu/j^j(Q)(i;) and 
p-p>v. Let qu, qv £ Q he the nodes of Q with which u and v were tagged, respectively. It follows 
that there is a fiow path from p to qu in the watershed overlay of Wjj^(m) and 'WR^J(Q){qu), so 
p G Wu (qu)- Analogously, p G Wy (qv)- Since {u,v) S X, we have qu 7^ qv, so any point p that 
is output by the algorithm lies in the intersection of the potential watersheds of two different 
nodes from Q. 

Next, we will argue that if p lies in the intersection of the potential watersheds of two 
different nodes from Q, then the algorithm will output p. Let q £ Q he the node with which 
p is tagged (hence, p G (q)), and let q' £ Q,q' ^ q he another node from Q such that 
p G Wu {q'). Consider a fiow path vr from p to q' in R\j{q'), and let (r, r') be the edge on vr 
such that r is tagged with a node other than q' and all nodes of TT[r',q'] are tagged with q' . 
Note that (r, r') must exist because all nodes of vr lie in Wy (Q) and have received a tag, p is 
tagged with another node than q' , and, since none of the nodes of Q lie in each other's potential 
watersheds, q' is tagged with itself. Therefore {r,r') exists, and (r, r') G X. Moreover, we have 
elev ji^(^Q-){r') = elev ji^(^qi-^{r'). Therefore r was put in the priority queue with the minimum 
elevation such that there is a realization R with elevf;{r') = clev j^^jf^qt-^ir') in which water fiows 
on the edge from r to r'. By induction on the nodes of tt from r back to p, it follows that p 
must eventually be extracted from the priority queue and output. 

This completes the proof of the lemma. □ 

5.5 The fuzzy watershed decomposition 

In this section we further characterize the structure of imprecise terrains by considering the 
ridge lines that delineate the "main" watersheds. In fact, the fuzzy watershed boundaries 
(Definition |8]) of the imprecise minima (Definition [5]) possess a well-behaved ridge structure if 
the terrain is regular. Consider the following definition of an "imprecise" ridge. 

Definition 9 Let Si, ..Sk be the imprecise minima of an imprecise terrain. We call the union 
of the pairwise intersection of the potential watersheds of imprecise minima the fuzzy ridge of 
the terrain. 

Let S he an imprecise minimum of a regular imprecise terrain. The next lemma testifies 
that the persistent watershed of any proxy g of 5 is equal to the intersection of the persistent 
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watersheds of all possible subsets of S. Therefore, we think of W(T^(q) as the actual minimal 
watershed of S, or the minimum associated with S. By Observation [T| the potential watersheds 
of all subsets of S are equal. Consequently, we think of the fuzzy watershed boundary of q as 
the fuzzy watershed boundary of S. 

Lemma 17 Let S be an imprecise minimum on a regular terrain, and let x be any proxy of S. 
Then r\Qcs^<^iQ)=^<^i^)- 

Proof: We want to argue about the intersection of the persistent watersheds of all subsets of S. 
Consider the complement C of this set, 

C:= I n I = U (Wn(Q)r= U w\«((Wu(Q)r). 

\QCS J QCS QCS 

By Observation [T] we have Wy (Q) = Wy (x) = (S) for any Q C S, so we have C = 
UqcsWI^*^ ((Wu iS)y). Now, the nodes contained in C can be characterized as follows. For any 
node p G C, there must be a realization, in which there is a subset S" C 5 and a node q outside 
Wu {S), such that there is a flow path from p to q that does not contain any node of S'. The 
given node {x} always serves as such a subset S' that is being "avoided", since x is a proxy and, 
by Definition [7| it is impossible for water that reaches x to continue to flow to a node outside 
of'Wu(5'). Therefore, 

((Wu {S)r) = Wu ((Wu {S)r) = W\J ((Wu {x)r) = C. 
The claim now follows from the definition of persistent watersheds. □ 

We can now further characterize the fuzzy ridge for regular terrains. The following lemma 
implies that on a regular terrain, the fuzzy ridge is equal to the union of the uncertainty areas 
of the fuzzy watershed boundaries of any representative set of proxies of the imprecise minima, 
see Corollary [1} 

Lemma 18 Let Si, ..Sk be the imprecise minima of a regular imprecise terrain and let qi, ..qk 
be associated proxies. For any 1 < i < k, we have that 'Wir]{qi) = ({Jj^i^u {<lj] 

Proof: Since qi is a proxy, we have that, 

{^MT = ^)^' ((Wu {q^)r) = Wu ((Wu {qm = [j Y^yj {p) . 

pe(Wufe))" 

Now, for a node p £ (Wu (^i))^, consider a minimum S that is reached by a flow path from 
p in . By Definition |6| we have that S is an imprecise minimum, and since p G (Wu {qi)y = 
(Wu {Si)y, we have S ^ Si. As such, S must be equal to some Sj for j ^ i. Furthermore, by 
Observation [1] we have W^-(S'j) = W^-(gj), therefore p G ^R-{(lj)- Now, Lemma 13 implies 
that Wu (p) C Wu {qj). It follows that (Wf^(gi))^ C [j^^. Wu {qj) . 

Since we also have that qj G (Wu {qi)Y for any j / i, we also get (Wpi(gj))'^ 5 Uj^j "^u {qj) > 
which implies the equality. □ 
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Corollary 1 Lemma 18 implies that, given qi,...,qig, a representative set of proxies for the 



imprecise minima of a regular imprecise terrain, it holds that 

U(Wu(gO\Wn(g,)) = U I UWu(g 

i i \ 

= U (^u(%)n|J Wu(g, 
i j^i 

By Observation^ this is equal to the fuzzy ridge of this terrain as defined in Definition^ This 




relationship is illustrated in Figure 12 
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Potential Watersheds 
Figure 12: Illustration to the fuzzy ridge on a regular terrain. 
Combining this with Lemma [9] and Lemma 16 we obtain: 



Theorem 6 We can compute the fuzzy ridge of a regular imprecise terrain in O(nlogn) time, 
where n is the number of edges of the imprecise terrain. 

Note that Definition [9] can also be applied to non-regular terrains, since it is solely based 



on the potential watersheds of the imprecise minima. We can use the algorithm of Section 5.2 



to compute proxies for these minima, and then use the algorithm of Lemma 16 to compute a 
fuzzy ridge between the watersheds of these proxies efficiently for non-regular terrains. However, 
note that the result may not be exactly the same as the fuzzy ridge according to Definition [9| 
because on a non-regular terrain, the potential watersheds of the proxies may be smaller than 
the potential watersheds of the imprecise minima. 



6 Conclusions 

In this paper we studied flow computations on imprecise terrains under two general models of 
water flow. For the surface model, where flow paths are traced across the surface of an imprecise 
polyhedral terrain, we showed NP-hardness for deciding whether water can flow between two 
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points. For the network model, where flow paths are traced along the edges of an imprecise 
graph, we gave efficient algorithms to compute potential (maximal) and persistent (minimal) 
watersheds and potential downstream areas. Our algorithms also work for sets of nodes and 
can therefore be applied to reason about watersheds of areas, such as lakes and river beds. 

In order to enable several extensions to these results in the network model, we introduced a 
certain class of imprecise terrains, which we call regular. We first defined when a set of vertices 
on an imprecise terrain can be considered a 'stable' imprecise minimum. We then described how 
to turn a non-regular terrain into a regular terrain using an algorithm by Gray et al. [11] and 
showed that this regularization algorithm preserves these imprecise minima. Interestingly, this 
algorithm also minimizes the number of minima of the terrain, while respecting the elevation 
bounds, as shown in |llj . 

We showed that persistent watersheds are nested on regular terrains and that these terrains 
have a fuzzy ridge structure which delineates the persistent watersheds of these stable minima. 
We gave an algorithm to compute this structure in O(nlogn) time, where n is the number of 
edges of the terrain. The correspondence between the imprecise minima of the regular and the 
non-regular terrain suggests that this fuzzy watershed decomposition on the regular terrain also 
allows us to reason about the structure of the watersheds on the original non-regular terrain. 
We think that, even though our work is motivated by geographical applications, the results will 
be useful in other application areas where watersheds are being computed, for instance in image 
segmentation [23J. 

There are many open problems for further research. Clearly, the contrast between the results 
in the surface model vs. the results in the network model leaves room for further questions, e.g., 
can we develop a model to measure the quality of approximations of water flow in the surface 
model, and how does it relate to the network model? Other flow models have been proposed in 
the GIS literature, e.g. D-c«, in which the incoming water at a vertex is distributed among the 
outgoing descent edges according to steepness. These models can be seen as modified network 
models which approximate the steepest descent direction more truthfully. In order to apply the 
techniques we developed for watersheds, we first need to formalize to which extent a node is 
part of a watershed in these models. 

Acknowledgments. We thank Chris Gray for many interesting and useful discussions on the topic 
of this paper. 

References 

[1] M. d. Berg, P. Bose, K. Dobrint, M. J. v. Kreveld, M. H. Overmars, M. d. Groot, T. Roos, 
J. Snoeyink, and S. Yu. The complexity of rivers in triangulated terrains. In Proceedings of the 8th 
Canadian Conference on Computational Ceometry, pages 325-330, 1996. 

[2] M. Borga, E. Gaume, J. Creutin, and L. Marchi. Surveying flash floods: gauging the ungauged 
extremes. Hydrological Processes, 22:3883-3885, 2008. 

[3] W. Buytaert, D. Reusser, S. Krause, and J. Renaud. Why can't we do better than Topmodel? 
Hydrological Processes, 22:4175-4179, 2008. 

[4] T. H. Gormen, G. E. Leiserson, R. L. Rivest, and G. Stein. Introduction to Algorithms. The MIT 
Press and McGraw-Hill Book Gompany, third edition, 2009. 

[5] W. Graddock, E. Kirby, N. Harkins, H. Zhang, X. Shi, and J. Liu. Rapid fluvial incision along the 
Yellow River during headward basin integration. Nature Geoscience, 3:209-213, 2010. 

[6] A. Danner, T. M0lhave, K. Yi, P. K. Agarwal, L. Arge, and H. Mitasova. TerraStream: from 
elevation data to watershed hierarchies. In Proc. 15th ACM Int. Symp. on Geographic Information 
Systems (ACM-GIS 2007), pages 212-219, 2007. 



31 



[7] M. de Berg, O. Cheong, H. Havcrkort, J.-G. Lim, and L. Toma. The complexity of flow on fat 
terrains and its i/o-efRcient computation. Computational Geometry, 43(4):331 - 356, 2010. Special 
Issue: 10th Workshop on Algorithms and Data Structures (WADS 2007). 

[8] M. de Berg, O. Cheong, H. Haverkort, J.-G. Lim, and L. Toma. The complexity of flow on fat 
terrains and its I/O-efRcient computation. Computational Geometry, 43(4): 33 1-356, 2010. 

[9] M. de Berg, H. Haverkort, and C. Tsirogiannis. Implicit flow routing on terrains with applications to 
surface networks and drainage structures. In Proc. 22nd ACM-SIAM Symp. on Discrete Algorithms 
(SODA), pages 285-296, 2011. 

[10] P. F. Fisher and N. J. Tate. Causes and consequences of error in digital elevation models. Progress 
in Physical Geography, 30(4):467-489, 2006. 

[11] C. Gray, F. Kammer, M. Lofflcr, and R. I. Silveira. Removing local extrema from imprecise terrains. 

CoRR, abs/1002.2580, 2010. 

[12] C. Gray, M. LofHer, and R. I. Silveira. Smoothing imprecise 1.5D terrains. Int. J. Comput. Geometry 
AppL, 20(4):381-414, 2010. 

[13] H. Haverkort and C. Tsirogiannis. Flow on noisy terrains: An experimental evaluation. In Proc. 
19th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems 
(ACM CIS), 2011. To appear. 

[14] F. Hebeler and R. Purves. The influence of elevation uncertainty on derivation of topographic 
indices. Geomorphology, 111(1-2) :4-16, 2009. 

[15] M. Henzinger, P. Klein, S. Rao, and S. Subramanian. Faster shortest-path algorithms for planar 
graphs. J. Computer and System Sciences, 55(l):3-23, 1997. 

[16] Y. Kholondyrcv and W. Evans. Optimistic and pessimistic shortest paths on uncertain terrains. In 
Proc. 19th Canad. Conf. on Com.put. Geom., pages 197-200, 2007. 

[17] R. D. Koster, S. P. P. Mahanama, B. Livneh, D. P. Lettenmaier, and R. H. Reichle. Skill in 
streamflow forecasts derived from large-scale estimates of soil moisture and snow. Nature Geoscience, 
3(9):613-616, August 2010. 

[18] J. Lindsay and M. Evans. The influence of elevation error on the morphometries of channel networks 
extracted from DEMs and the implications for hydrological modelling. Hydrological Processes, 
22(11):1588-1603, 2008. 

[19] Y. Liu and J. Snoeyink. Flooding triangulated terrain. In Proc. 11th Int. Symp. Spatial Data 
Handling, pages 137-148, BerUn, 2005. 

[20] M. Mcallistcr and J. Snoeyink. Extracting consistent watersheds from digital river and elevation 
data. In In Proc. ASPRS/ACSM Annu. Conf, 1999. 

[21] A. Montanari. What do we mean by 'uncertainty'? The need for a consistent wording about 
uncertainty assessment in hydrology. Hydr. Proc, 21:841-845, 2006. 

[22] J. O'Callaghan and D. Mark. The extraction of drainage networks from digital elevation data. 
Computer vision, graphics, and image processing, 28(3):323-344, 1984. 

[23] D. L. Pham, C. Xu, and J. L. Prince. Current methods in medical image segmentation. Biomedical 

Engineering, 2:315 337, 2000. 

[24] B. Sivakumar. The more things change, the more they stay the same: the state of hydrologic 
modelhng. Hydrological Processes, 22:4333-4337, 2008. 

[25] D. Tarboton. A new method for the determination of flow directions and upslope areas in grid dig. 
elev. models. Water Resources Research, 33(2):309-319, 1997. 

[26] D. Tetzlaff, J. McDonnell, S. Uhlcnbrook, K. McGuire, P. Bogaart, F. Naef, A. Baird, S. Dunn, 
and C. Soulsby. Conceptualizing catchment processes: simply too complex? Hydrological Processes, 
22:1727-1730, 2008. 



32 



[27] J. A. Vnigt, C. G. H. Diks, H. V. Gupta, W. Boutcn, and J. M. Vcrstratcn. Improved treatment 
of uncertainty in hydrologic modeling: Combining the strengths of global optimization and data 
assimilation. Water Resources Research, 41, 2005. 

[28] S. P. Wechsler. Uncertainties associated with digital elevation models for hydrologic applications: 
a review. Hydrology and EaHh System Sc., 11(4):1481-1500, 2007. 



33 



A Computing potential watersheds in linear time 



Theorem 3 The canonical realization of the potential watershed of a set of cells Q in an im- 
precise grid terrain of n cells can be computed in 0{n) time. 



Proof: The computation of potential watersheds in Section 4.1 has much in common with 
computing single-source shortest paths. In both cases, the goal is to compute a label 5{v) for 
each node v: in the case of potential watersheds it is the lowest elevation such that a flow path 
to a given destination q exists; in the case of shortest paths it is the distance from the given 
source q. During the computation, we maintain tentative labels d[v] for each node v which are 
upper bounds on the labels to be computed. (The tentative label of a node that has not been 
discovered yet would be c«.) The computations consist of a sequence of edge relaxations: when 
relaxing a directed edge {u,v), we try to improve (that is, lower) d[v] based on the current 
value of d[u], which is an upper bound on S{u). Both problems share some crucial properties: 
for every node v that can be reached, there is a "shortest" path ■k{v) = uo,ui, ...,Uk where 
uq = q and Uk = v, the correct labels S{uo), 5{ui), 6{uk) form a non-decreasing sequence, and 
when the edges on this path are relaxed in order from {uo,ui) to {uk-i,Uk), the relaxation of 
{ui-i,Ui) will correctly set d[ui] equal to 6{ui). All that is necessary for the computations to 
compute all labels, is that the sequence p of relaxations performed by the algorithm contains 
7r{v) as a subsequence, for each v. Note that the edges of Tr{v) do not need to be consecutive in 
p: the labels along tt{v) are computed correctly even if the relaxations of 7r{v) are interleaved 
with relaxations of other edges, or even with out-of-order relaxations of edges of Tr{v). 

There are several algorithms to find a sequence of relaxations p in the above setting, such 
that for every node v, the sequence p contains the relaxations of a shortest path tt{v) as a 
subsequence. These algorithms are usually known as algorithms to compute (single-source) 
shortest paths, but they can also be applied directly to the more general setting described 
above. Dijkstra's algorithm finds a sequence of relaxations that is optimal in the sense that it 
relaxes each edge only once. However, to achieve this, the algorithm needs 0(n) operations on 
a priority queue of size Q{n) in the worst case, where n is the number of nodes and edges in 
the graph P]. 

An alternative is the algorithm of Henzinger et al. pL5jj . This algorithm uses a hierarchy of 
priority queues. Most priority queue operations in this algorithm are on small priority queues. 
The algorithm needs more relaxations than Dijkstra's algorithm, but still not more than 0(n). 
Provided the relaxations take constant time each, the whole algorithm runs in 0(n) time. 
However, the algorithm by Henzinger et al. only works if a recursive decomposition of the graph 
is provided that satisfies certain properties. Fortunately such decompositions can be found in 
0{n) time for planar graphs, and also for certain other types of graphs. In particular, it is easy 
to construct such a decomposition for a graph that represents a grid terrain model, even in the 
model where each cell can drain to one or more of its eight neighbors, for which the adjacency 
graph is non-planar. Let ri < r2 < ... be a sequence of powers of four. Now we can easily 
make a decomposition of the graph into square regions of x .yri nodes; we group these 
together into regions of x regions, etc., generally grouping regions of x y^ nodes 
into regions of y^^i+i x y^ri+i nodes (some regions at the boundary of the whole input grid may 
be slightly smaller). On each level i, the regions have size 0(rj) and each region has @{^/ri) 
nodes on its boundary, thus each level forms a so-called rj-division. We choose the region sizes 
such that they satisfy Equation (19) from Henzinger et al. 

With this decomposition, the structure of the single-source shortest paths algorithm from 
Henzinger et al. can also be applied to the computation of potential watersheds on grid terrains. 
For grid terrains, dmax = 0(1), and thus, the computation of the slope diagrams and the 0{n) 
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relaxation steps from the "shortest-paths" algorithm take only 0{n) time. Together with 0(n) 
time for priority queue operations, we get a total running time of 0(n). □ 

B Persistent watersheds with multiple connected components 

Lemma 19 There exists a regular terrain that contains a persistent watershed that consists of 
more than one connected component. 

Proof: Refer to Figure [TS} The figure shows five nodes with their elevation intervals. The 
edges {a,b), {b,d) and (c, d) have length 1. The edge (d, e) has length 1.6. From a and e, very 
steep edges lead to nodes downwards not shown in the figure. The potential watershed Wy (e) 
of e is {c, d, e}. The node d is not in the persistent watershed of e: if d has elevation more than 
6|, the flow path from d will lead to b, outside Wy (e). In that case c is a local minimum inside 
Wy (e). Whenever c is not a local minimum, the elevation of d must be less than 4, and the 
flow path from c will lead to d and on to e. Thus c is in the persistent watershed Wfr|(e) of e, 
but d is not, so we have WfT|(e) = {c, e}. □ 

c [4,4] 
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Figure 13: Example of disconnected persistent watershed on a regular terrain. 
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